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SECTION  I 


INTRODUCTION  AND  SUMMARY 

Computation  of  high- resolution  wavenumber  spectra  includes, 

as  an  intermediate  step,  the  solution  of  a  set  of  Hermi.ian  equations.  This 

set  of  equat;  ms  has  the  form  of  a  least-squares  multichannel  frequency- 

domain  filter  design  equation  which  can  be  expressed  in  vector-matrix 
notation  as 


H  a  =  b 


where 


(i-D 


H  =  n  x  n  nonsingular  Hermitian  power  spectral  matrix 

b  =  nx  lknowncomplexcolumnvectorofoutputpowerspectra 
n  x  1  “"known  complex  column  vector  of  filter  weights 


This  report  investigates  three  techniques  of  solving  for  the  unknown  vector 

a:  the  method  of  conjugate  gradients,  steepest-descent  method,  and  exact- 

inverse  method.  The  object  is  to  determine  the  accuracy  and  computational 
complexity  of  each  technique. 


Previously,  problems  of  this  type  were  solved  using  a  direct 

numerical  matrix  technique  such  as  the  Gaussian  elimination  method4  or  the 
square -root  method. 


Fox,  L  1965:  An  Introduction  to  Numerical  Linear  Algebra 
Oxford  University  Press,  Nr-w  York,  p.  205-213. 


Faddeeva,  V.N.  ,  1959:  Computational  Methods 

Dover  Publications,  Inc.,  New  York. 


of  Linear  Algeb 


ra. 
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To  obtain  a  solution,  the  Gaussian  elimination  method  requires 
only  that  the  matrix  be  nonsingular,  while  the  square- roots  method  requires 
the  matrix  to  be  nonsingular  and  Hermitian.  It  was  apparent  that  neither  a 
knowledge  of  the  form  of  the  H  matrix  nor  an  estimate  of  the  solution  vector  a 
would  significantly  simplify  the  computation  involved  with  either  of  these  tech¬ 
niques.  Since  both  the  general  form  of  the  H  matrix  and  an  estimate  of  the 
solution  vector  were  available,  an  investigation  of  two  iterative  techniques 
(steepest-descent  and  conjugate-gradients)  which  appeared  to  benefit  compu¬ 
tationally  from  this  knowledge  was  initiated.  Found  during  this  investigation 
was  a  theorem  from  linear  algebra  which  analytically  expresses  the  inverse 
of  the  particular  H  matrix  used  to  generate  high- resolution  f-k  filter  sets. 
Application  of  this  analytical  inverse  greatly  reduced  the  computation  needed 
to  generate  high- resolution  f-k  filter  sets. 

The  high- re  solution  wavenumber  filter  sets  are  designed 
using  a  power  spectral  matrix  H  of  the  form 

h  =  [«i  +  ;;*]  (i.2) 


where 

I  =  n  x  n  identity  matrix 

x  =  n  x  1  column  vector  of  the  channel  transforms 
— # 

x  =  n  x  1  row  vector  of  the  conjugate  channel  transforms 
Using  the  exact  inverse  equation,  the  filter  vector  a  can  be  computed  easily  as 


(1-3) 


1-2 


•ol«nc«  ••rvio 


llvli 


This  computation  requires  much  fewer  calculations  than  a  single  cycle  of  the 
iterative  techniques  studied  and  is  orders  of  magnitude  simpler  than  the 
Gaussian  elimination  technique. 

Further  application  of  the  exact  inverse  equation  arises  when 

P  spectral  matrix  is  generated  using  an  exponentially  weighted  series 
of  transform  vectors: 


N 


H 


-  E 

j=i 


N-j - * 

a  J  x.  x 
J  J 


0  <  a  £  1 


(1-4) 


The  inverse  of  H  can  be  generated  iteratively  by  writing 


Hi+1  =  a  +  x  x* 
J+1  J  J+l  J+l 


d-5) 


The  exact  inverse  equation  is 


H 


-1 


j+l 


a 


H 


-1 


1 


—  *  .  l  — 

a  +  X.x]  H.  X.  , 
J+1  J  j+l 


TT-1-  -*  .  1 

H-  X  X  H. 

J  J+l  J+l  J 


(1-6) 


This  formulation  for  the  inverse  of  a  power  spectral  matrix  is  especially 
useful  when  the  transform  vectors  X.  are  available  sequentially  as  they  would 
be  in  a  practical  situation.  The  inverse  of  the  power  spectral  matrix  can  be 
URdated  as  each  transform  vector  becomes  available!  and  from  this  inverse 
matrix,  many  frequency-domain  filter  sets  can  be  designed  at  each  iteration. 
The  ability  to  track  a  spectral  matrix  adaptively  facilitates  the  coherent  fre¬ 
quency-domain  processing  of  small  arrays  with  nonstationary  noise  fields. 
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SECTION  II 


MULTICHANNEL  FREQUENCY  FILTER  DESIGN 

This  itport  investigates  three  iterative  techniques  which  offer 

computational  savings  in  the  design  of  multichannel  frequency-domain  filter 

sets.  This  study  was  initiated  to  reduce  the  computational  complexity  involved 

— ♦ 

in  generating  high- resolution  f-k  spect.-a,  thereby  enhancing  the  technique's 

capability  of  a  real-time  detection  and  location.  The  following  discussion 
briefly  considers  the  general  multichannel  frequency-domain  design  problem 
then  proceeds  to  the  special  problem  of  high- resolution  f-k  filter  design. 

The  general  •nultichannel  filter- design  equation  can  be  written 
in  vector-matrix  notation  as 


X1  X1 


x  x. 
n  1 


x,  x 
1  n 


* 

x  x 
n  n 


n 


•* 

X1S 


■f 

x  S 
n 


=  b 


(2-1) 


$ 

where  x.  x.  is  the  estimated  crosspower  (for  i  £  j)  or  the  estimated  autopower 

1  ^  th  th  * 

(for  i  =  j)  between  the  i  and  the  j  channels  and  where  x.  S  is  the  cross- 

J|>h 

power  estimate  between  the  desired  signal  output  and  the  j  data  channel. 
These  autopower  and  crosspower  spectra  usually  are  estimated  from  a  sample 
of  multichannel  data  using  one  of  two  general  approaches. 
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One  method  is  generating  the  autocorrelations  and  cross  cor¬ 
relations  from  the  time-domain  data  segment  and  transforming  these  cor¬ 
relations  to  obtain  the  desired  power  spectra.  The  other  method  involves 
transforming  segments  of  the  data  sample  and  stacking  and/or  smoothing  these 
transforms  to  obtain  the  desired  power  spectra.  The  second  method  is  perti¬ 
nent  to  this  report  as  it  is  computationally  more  efficient  than  the  first  method. 


vector 


The  transform  operation  itself  is  not  discussed,  but  the  data 


(2-2) 


j.s  assumed  to  be  the  colateral  discrete  transforms  of  all  the  data  channels  of 

the  -f/t*1  data  segment  at  a  particular  frequency.  The  power  spectral  matrix 

at  this  frequency  then  can  be  estimated  as  the  weighted  sum  of  a  priori  signal 

- * 

power  spectral  information  S  and  the  transform  crossproducts  x  x  so  that 


m 


=  H  =  i  S  +  £ 


t=l 


—  —  * 


(2-3) 


* 

x  x 
n  n 
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Using  this  estimate  for  the  power  spectral  matrix,  the  filter  design  equation 


can  be  written  as 


m  _  _  "I 

«s  +  E  nxi  xi\  a  =  b 

i  -  i  J 


(2-4) 


Most  iterative  solution  techniques  compute  a  residual  vector 


r^  from  an  initial  guess  v^  of  the  solution  vector  a: 


rl  =  b  •  H  vo 


(2-5) 


In  general,  the  computation  of  this  residual  vector  requires  n  complex 
multiply  and  add  operations,  not  including  the  computations  needed  to  form 
the  matrix  H-  In  the  special  case,  using  Equation  2-3  to  form  the  matrix  H, 
the  residual  vector  can  be  computed  directly  without  explicitly  forming  the 
H  matrix.  This  is  done  by  substituting  Equation  2-3  into  Equation  2-5  and 
using  a  different  order  of  computation,  as  shown  in  the  following  equation. 


[m  _  -i 

l8t£B*=*:‘]v 


(2-6) 


5  S  v0  +  E 


ai  Xl 


(**  ’«)  - 


To  form  the  H  matrix  using  Equation  2-3,  mn  complex  multi¬ 
ply  and  add  operations  are  required.  Equation  2-6  can  compute  the  residual 
2 

using  n  +  3mn  complex  multiply  and  add  operations  by  first  forming  the 

—  — 

vector  dot  product  x^  vQ  ,  then  performing  the  scalar  multiplication,  and 

2 

summing  over  the  resulting  vectors.  This  compares  to  (m  +  l)n  complex 
multiply  and  add  operations  needed  to  form  H  separately  and  then  compute  r. 
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Further  Dmputational  savings  can  be  obtained  if  the  signal- 

2 

model  matrix  S,  which  contributes  the  n  term,  can  be  expressed  as  a  diago¬ 
nal  matrix  (uncorrelated  noise) 


0 


S  = 


0 


0 


(2-7) 


or  the  sum  of  vector  products  (several  plane  waves) 


S 


'k  k  k 


(2-8) 


For  diagonal  matrix  S,  (3m  +  1)  n  complex  multiply  and  add  operations  are 
required  to  compute  r.  If  S  were  defined  by  Equation  2-8,  3  (m  +  p)  n 

multiply  and  add  operations  would  be  needed.  Since  using  the  Gaussian 
elimination  method  to  solve  Equation  2-1  requires  at  least  2n  operations, 
significant  savings  in  computation  will  result  from  using  iterative  solution 
techniques,  provided  only  a  few  iterations  are  required. 

The  special  problem  of  designing  high- resolution  f-k  filter 
sets  uses  a  simple  form  of  the  power  spectral  matrix  where 

H  =  [n  +  xx"]  (2-9) 
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This  form  results  from  the  transient  nature  of  the  detected  signals,  but 
these  likely  will  not  be  present  in  several  different  transform  segments 
of  the  data  sample.  The  signal  model  is  assumed  to  be  equalized  and 

uncorrelated  from  channel  to  channel;  therefore,  S  =  i  (the  identity  matrix), 
and 


1 

0 


b  = 


(2-10) 


A  theorem  from  linear  algebra,  which  will  be  called  the 
exact  inverse  equation  (Appendix  B),  reduces  the  computation  of  high- 
res  olutic  n  f-k  filter  sets  to  approximately  3n  complex  multiply  and  add 
operations  (where  n  is  the  number  of  channels).  The  exact  inverse 
equation  can  be  expressed  generally  as 


where 


[A  +  uv*] 


-1 


=  A 


-1--*  - 1 

A  u  v  A 

-*  -1- 
1  +  v  A  u 


A  =  n  x  n  nonsingular  complex  matrix 
u  =  n  x  1  column  vector 

v  =  n  x  1  column  vector 
v  =  conjugate  transpose  of  v 


[a  +  u  ;*]i 


is  nonsingular 


(2-11) 
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Applying  this  theorem  to  the  high- resolution  f-k  filter 


design  yields 


a 


- # 


x  x 


— *  - 
+  X  X 


b 


(2-12) 


Pointing  out  the  simplicity  of  the  calculations,  this  vector- matrix  equation 
can  be  rewritten  as 


a 


1 


a. 

J 


(2-13) 


—  *  — 

The  vector  product  x  x  is  real  and  needs  to  be  computed 

only  once  for  each  transform  segment.  Remaining  calculations  can  be 

accomplished  with  the  equivalent  of  2n  more  complex  multiply  and  adds 

(CMPA's).  The  required  computations  are  discussed  in  Section  V.  The 

total  computation  necessary  for  the  exact  inverse  method  is  roughly  3n 

complex  multiply  and  add  operations  compared  with  approximately  12n  CMPA's 

for  one  iteration  of  the  steepest-descent  method  and  24n  CMPA's  for  one 

iteration  of  the  finite  procedure.  Thus,  the  exact  inverse  technique  is  the 

most  straightforward  and  efficient  procedure  to  design  high- resolution 
— ♦ 

f-k  filter  sets.  It  is  possible  to  manipulate  the  exact  inverse  equation  to 
update  the  inverse  of  a  power  spectral  matrix;  this  area  is  covered  in  the 
section  on  adaptive  update  methods. 
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SECTION  III 

SUMMARY  OF  GRADIENT  METHODS 

Based  on  the  gradient  of  a  particular  quadratic  form,  two 
iterative  design  procedures  —  steepest-descent  and  conjugate- gradient 
methods  -  were  investigated  during  this  study.  While  these  two  approaches 
are  not  the  most  efficient  means  of  obtaining  high- resolution  f-k  filter  sets, 
they  are  applicable  to  the  more  general  filter-design  problems.  The  details 
of  these  methods  are  dtscribed  in  Appendix  A. 

A.  STEEPEST-DESCENT  METHOD 

The  numerical  algorithm  generated  using  the  steepest- 
descent  approach  is 

"k+l=\+Pk7k  (3-1) 

(3-2) 

(3-3) 


A  simpler  procedure  results  when  the  form  of  the  H  matrix 
is 

m 

H  =  * 1  +  £  <3-4> 
t= i 

The  matrix  operations  become  vector  operations,  and  the  algorithm  is 


=  b  -  H  v 


Pi.  = 


—  *  _ 
ri  ri 

k  k 

—  *  _ 
r,  H 
k  k 


Vk+1  =  Vk  +  Pkrk 


(3-1) 


m 


rk  =  b  "  *  Vk  '  T,  Wl  *1 
1=  1 


(x^  vk) 


(3-5) 
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pk  = 


—  *  — 
rk  rk 


— *  _ 


5  rk  rk  +  £  » 


(3-6) 


l  Zl  Zl 


-#  - 
Xl  rk 


(3-7) 


While  the  notation  of  Equations  3-5,  3-6,  and  3-7  is  more 

complicated  than  that  of  Equations  3-2  and  3-3,  the  amount  of  computation. 

needed  has  been  reduced.  The  reduction  in  computation  is  obtained  by  per- 

—  $  —  —  # 

forming  the  vector-dot-product  operations  x.  v  and  x.  r  first,  leaving 
only  scalar  operations  to  be  performed.  An  analysis  of  the  amount  of  com¬ 
putation  required  is  presented  in  a  later  section  of  this  report. 

B.  FINITE  PROCEDURE  (CONJUGATE- GRADIENT  METHOD) 

The  computation  cycle  evolving  this  procedure  is 


=  v  +  p  w, 
k  k  k 


(3-1) 


=  b  -  H  v. 


—  *  — 
r  r 
k  k 

-*  - 

r.  H  r, 

k  k 


(3-2) 


(3-3) 


-#  — 
rk  rk 

*k-l  =  - 

rk-l  rk-l 


(3-8) 


Wk  =  rk  +  Vl  Wk-1 


(3-9) 


The  first  cycle  is  computed  using  =  r^. 


Ill- 2 


aolsnct  ssrvleos  division 


(3-4) 

(3-1) 

(3-5) 

(3-10) 

(3-7) 

(3-11) 

(3-12) 


qk-  1 


(3-9) 


For  the  first  iteration,  wrt  =  r.  and  q,  ,  is  not  computed. 

0  U  k-1 

Notation  again  becomes  more  cumbersome,  but  the  computation 
needed  is  much  less  than  that  required  to  form  the  autopower- crosspower 
matrix  and  then  solve  for  the  desired  filter  set.  Computational  savings  are 
obtained  by  carrying  out  the  same  computations  in  a  different  order.  A 
detailed  analysis  of  the  numerical  operations  involved  is  presented  in  a  later 
section  of  this  report. 
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SECTION  IV 

ADAPTIVE  UPDATE  METHODS 


In  many  practical  applications,  all  individual  vectors  x. 
which  make  up  the  estimated  covariance  matrix  H  are  not  simultaneously 
available.  The  various  samples  must  be  stored  before  they  can  be  used 
in  the  filter  design.  One  method  stores  all  the  needed  vectors  in  a  com¬ 
puter  memory  and  uses  one  previously  mentioned  algorithm  to  compute 
a  new  filter  set  each  time  a  new  vector  is  received.  Another  method 
iteratively  updates  either  the  filter  set  or  the  covariance-matrix  inverse 
each  time  a  new  vector  is  received.  Iterative  updating  reduces  the  storage 
and  computation  required  but  is  restricted  to  using  a  class  of  memory 
functions.  This  class  includes  the  exponential-weighting  function  which  is 
discussed  in  this  section. 

Two  of  the  techniques  previously  discussed  can  be  formulated 
as  iterative -updating  techniques  which  weight  previous  transform  vectors 
with  an  exponentially  decreasing  function.  Exponential  weighting  can  be 
obtained  by  computing  a  weighted  sum  of  the  previous  covariance  estimate 
and  the  current  transform-vector  product;  i.  e.  , 

H  =  [a  H  +  x  x  '  1  Osa^l  (4-1 

n  |_  n- 1  n  n  J  ' 


etc.  Thus,  a  vector  n  samples  in  the  past  is  weighted  by  a".  If  0  <  a  <  1, 
the  current  estimate  of  the  autopower- crosspower  matrix  weights  the 
current  transform  vector  with  weights  equal  to  unity  and  previous  transform 
vectors  with  weights  of  less  than  unity.  These  weights  (Equation  4-2)  become 
exponentially  smaller  for  the  more  distant  samples.  The  past  transform 
vectors  can  be  weighted  in  the  desired  fashion  by  the  choice  of  cl,  where 
a  =  exp (-3)  for  0  >  0. 

A.  EXACT  INVERSE  MATRIX  UPDATE  METHOD 


One  updating  method  applies  the  exact  inverse  equation  directly 
to  Equation  4-1.  The  inverse  of  the  updated  covariance  matrix  is  then 


-1 

1 

H-1  x  x*H 
n-1  n  n  n-1 

n 

a 

—  *  T  - 1  — 

a  +  x  H  x 

n  n- 1  n 

(4-3) 


The  desired  filter  weights  can  be  obtained  from  the  matrix  multiplication 


f  = 


(4-4) 


In  this  case,  the  previously  computed  covariance  matrix  is  considered  to 
be  the  noise  model.  This  method  requires  computing  and  storage  of  an 
entire  matrix  at  each  iteration  as  well  as  a  matrix- vector  multiplication  for 
each  filter  set  needed.  When  only  a  few  filter  sets  are  needed,  the  steepest- 
descent  technique  is  computationally  more  efficient. 

B.  STEEPEST-DESCENT  METHOD 

A  particularly  simple  scheme  is  afforded  by  extending  the 
steepest-descent  method.  The  new  filter  set  is  formed  from  the  previous 
filter  set  by  computing 
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n+1  n  n+1  n+1 


where  the  residual  vector  r  ,  is 

n+1 


b  “  H  xi  f 
n+1  n 


The  matrix 


_  _ # 

a  H  +  x  x 
n  n+1  n+1 


This  allows  the  residual  vector  to  be  expressed  as 


rn+l  =  b 


'  [ttHn 


+  Xn+1  Xn+1 


] 


=  (l-a)b+afb-H  f  1-x  x*  .  f 
L  nnJ  n+1  n+1  n 


The  quantity  b  -  H  f  is  the  residual  vector  r  from  the  previous 

n  n  n 

iteration.  The  new  residual  vector  is  then 


r  ,  1  =  (l  -  a)  b  +  a  r  -xx  f 

n+l  n  n+1  n+1  n 


This  new  residual  vector,  used  to  calculate  the  new  filter  weights,  i; 
stored  to  be  used  in  the  next  residual  calculation. 
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The  problem  of  computing  the  iteration  factor  p  is  en- 

n+1 

countered  now.  Since  the  matrix  H  exists  implicitly  at  each  iteration  but 
is  not  directly  available  for  computation,  the  iteration  factor  cannot  be 
computed  directly  from 

—  *  - 
r  r 
nH  n+1 

Pn+ 1  ~  *  — 

r  H  r 
n+1  n+1  n+1 


An  alternative  to  computing  the  iteration  factor  for  each  cycle  is  to  choose 
a  single  value  and  use  it  throughout  the  process.  This  single  value  must  be 
small  enough  to  guarantee  the  numerical  stability  of  the  process.  This 
requirement  will  be  satisfied  if 

0  <  p  5  - - - 


where  is  the  largest  eigenvalue  of  the  matrix  Since  the  matrix 

is  changing  continually  as  new  data  samples  are  added,  only  the  upper 
bound  on  the  maximum  eigenvalue  can  be  expressed.  A  matrix  of  the  form 
of  Equation  4-1  is  essentially  a  weighted  sum  of  vector  products 


H  =  Y''  aJ  x  .x 
n  Cmd  n-j  n-j 

j=0 


for  0  <  a  <  1 


The  largest  eigenvalue  \n  of  H  is  then 

max  n 


\  < 

max 


Fx*  x  1 

|_  max  maxi 
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where  x  is  the  largest  data  vector  encountered. 

IllaX 

choice  of  p  should  be  in  the  range 


This  implies  that  the 


0<ps 


Such  a  conservative  choice  of  the  iteration  factor  sacrifices 

rapid  convergence  for  numerical  stability.  This  approach  requires  a  priori 

knowledge  of  the  maximum  magnitude  of  the  transform  vector.  In  a  practical 

situation,  overestimating  the  value  of  x *  x  would  be  necessary  to 

max  max  7 

maintain  p  within  the  desired  bounds. 

Another  approach  to  determining  pis  to  compute  x  x  at  each  cycle 

—  sjc  — 

and  compare  it  with  the  largest  previous  valueofx  x.  Thenpcanbe  chosenwithin 
the  bounds  determined  by  the  previous  data  us  ed  in  the  process .  This  approach 
might  permit  the  use  of  an  iteration  factor  which  is  initially  much  larger  than  can 
be  used  when  the  factor  is  fixed  throughout  the  process. 

C.  COMMENTS 

Neither  of  the  preceding  adaptive  techniques  have  been  applied 
to  either  real  or  synthetic  data.  A  simplification  of  the  steepest-descent 
algorithm,  similar  to  the  technique  investigated  by  Bernard  Widrow,  has 
been  investigated  using  time-domain  data.  That  investigation  is  described 

sjcsje 

in  Advanced  Array  Research  Special  Report  No.  1. 


Widrow,  Bernard,  1966,  Adaptive  Filters  1:  Fundamentals:  Stanford 
University  Tech.  Rpt.  No.  6764-6,  Contracts  DA-01-021  AMC-90015  (Y) 
and  NOBsr-95038,  Dec. 

Texas  Instruments  Incorporated,  1967:  Adaptive  Filtering  of  Seismic 
Array  Data,  Advanced  Array  Research  Spec.  Rpt.  No.  1,  Contract 
F33657-67-C-0708-P001,  18  Sept. 
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SECTION  V 

REQUIRED  COMPUTATIONS 


A  major  consideration  in  choosing  a  numerical  method  to 
solve  a  particular  problem  is  the  number  of  computations  required  by  each 
method.  This  section  presents  a  summary  of  the  computations  required  by 
each  method  discussed  in  this  report.  For  comparison,  the  method  of 
Gaussian  elimination  is  also  included.  An  example  is  given  at  the  conclusion 
of  the  discussion  to  illustrate  the  numerical  efficiency  of  each  method. 

Each  method  is  evaluated  on  the  basis  of  solving  the  following 
equation  for  the  elements  a.  of  the  vector  a. 


A.  GAUSSIAN  ELIMINATION  METHOD 


By  performing  elementary  row  operations  on  the  H  matrix 
and  the  b  vector.  Equation  5-1  is  changed  to  an  equation  of  the  form 
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.  ...  ..«*,*?* 


The  value  of  is  then  d^.  The  remaining  elements  of  a  are  found  by 

successively  solving  for  a  . ,  a  ,,  ....  a,.  To  use  this  method,  the 

n-i  n-2  1 

covariance  matrix  H  must  be  computed  directly.  This  requires  2n  multi- 

2 

plies  and  2n  adds  to  calculate  diagonal  elements  and  2n  -  2n  multiplies  and 
2 

n  -  n  adds  to  calculate  the  off-diagonal  elements. 

Then  the  matrix  is  reduced  to  the  triangular  form  of 
3  2  3  2 

Equation  5-2  which  requires  2n  -  2n  multiplies,  2n  -  2n  adds,  and 
2 

2n  divides.  Solving  for  the  filter  weights  using  the  triangular  form 

2  2 

(Equation  5-2)  requires  an  additional  2n  -  2n  multiplies  and  2n  -  2n  adds. 

For  all  the  computations,  the  Gaussian  elimination  method, 
including  the  covariance  matrix  formation,  requires 

3  2 

•  2n  +  2n  -  2n  multiplies 

3  2 

•  2n  +  n  -  n  adds 

2 

•  2n  divides 

B.  STEEPEST-DESCENT  METHOD 

< 

i . 

Calculating  the  residual  vector  requires  lOn  multiplies  and 
8n  +  1  adds.  To  form  the  new  filter  weights,  2n  multiplies  and  2n  adds 
are  required.  If  the  algorithm  is  cycled  for  ITER  iterations,  calculations 
using  the  steepest-descent  method  require 

•  12n  (ITER)  multiplies 

•  (lOn  +  1)  (ITER)  adds 
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C.  FINITE  PROCEDURE 

Assuming  the  imaginary  part  of  p  is  0,  each  iteration  of 

1C 

the  finite  procedure  uses  lOn  multiplies  and  8n  +  1  adds  to  compute  the 
residual  vector.  Computing  the  iteration  constant  requires  lOn  +  3  multi¬ 
plies  and  8n  +  2  adds  for  the  first  iteration  and  12n  +  3  multiplies, 
lOn  +  2  adds,  and  1  divide  for  each  succeeding  iteration.  To  calculate 
the  new  filter  weights,  2n  multiplies  and  2n  adds  at  each  iteration  are 
required.  If  there  are  ITER  iterations,  the  total  computation  requires 

•  (24n  +  3)  (ITER)  -  2n  multiplies 

•  (20n  +  3)  (ITER)  -  2n  adds 

•  ITER  divides 

D.  EXACT  INVERSE  METHOD  FDR  SINGLE- CHANNEL  PREDICTION 

Computing  the  filter  weights  by  this  method  requires 

•  3  n  multiplies 

•  n  +  1  adds 

•  1  divide 

E.  COMPARISON 

As  an  example,  consider  the  design  of  a  high- resolution  f-k 
spectra  filter  set  where  a  first  approximation  of  the  filter  weights  is  avail¬ 
able.  Using  20  channels  of  data,  the  required  amount  of  computation  for 
each  method  is  as  follows: 

•  Gaussian  Elimination  —  16, 760  multiplies, 

16,380  adds,  and  800  divides 

•  Steepest  Descent  (2  iterations)  —  480  multiplies 
and  402  adds 
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•  Finite  (1  iteration)  —  443  multiplies, 
363  adds,  and  1  divide 


•  Exact  Inverse  —  40  multiplies,  21  adds, 
and  1  divide 

In  this  example,  the  exact  inverse  procedure  is  clearly  the 
most  efficient  computationally.  This  efficiency  is  due  to  the  particular 
formulation  of  the  equations  being  solved,  and  the  solution  of  a  more 
general  filter-design  problem  mi<vht  show  one  of  the  other  techniques  to 
be  more  effective.  The  selection  of  a  technique  depends  on  the  problem 
formulation  and  on  the  knowledge  of  the  form  of  the  solution. 


SECTION  VI 
CONCLUSIONS 


Analysis  of  the  methods  available  for  the  numerical  solutior. 
of  the  frequency-domain  multichannel  filter-design  problem  yields  two 
major  conclusions: 


•  The  exact  inverse  equation  is  by  far  the 
most  satisfactory  method  for  designing 
high -re solution  filter  sets  from  single 
transform  data  (a  rank-one  matrix  of 
data). 

•  The  exact  inverse  equation  can  be  used 

to  update  the  inverse  of  a  spectral  matrix 
for  adaptively  tracking  nonstationary 
noise  fields.  Such  information  would  be 
required  to  do  Baysian  location  in  a  cor¬ 
related  noise  field  {for  example,  at  the 
subarray  level). 


APPENDIX  A 
GRADIENT  METHODS 


APPENDIX  A 
GRADIENT  METHODS 


A  method  of  iteratively  solving  the  equation  H  a  =  b  for  an 
unknown  a  can  be  generated  by  minimizing  the  quadratic  form 


~  — *T  -  — #  T  —  — *  T  — 

Q  =  v  Hv-b  v-v  b 


(A-  1 ) 


__  ^  rp 

When  H  is  positive -definite  Hermitian,  the  quadratic  v  Hv  is  non- 
negative,  and  Q  takes  on  its  minimum  value  -  b  Hb  when  v  =  H~  Id. 

Thus,  the  vector  v  which  minimizes  Q  is  the  solution  vector  a.  The  iter¬ 
ative  procedure  makes  an  improved  estimate  v^  of  the  solution  vector  by 
starting  with  an  initial  estimate  v^  and  successively  computing  a  new  estimate 


v,  ,  i  =  v  +  p  w, 
k+1  k  rk  k 


for  p  real 
k 


(A-2) 


The  amount  Q  changes  from  the  k  to  the  k  +  1  iteration 


n  n  2  -*  T  TT  -  /— *T  —  -*  T  -\ 

Qk+i '  Qk  -  pk  w  Hw-Mwk  \  +  rk  H 


(A-3) 


where 


r  =  b  -  Hv 
k  k 


Minimizing  Q  requires  that  AQ  be  negative  after  each  iteration.  This  con¬ 
strains  pk  to  the  range 


°  <  pk  < 


—  *T  —  —  *  T  — 

wk  rk  *  rk  wk 

—  *  T  — 
w,  H  w, 

k  k 
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The  largest  change  in  AQ  is  obtained  when 


P 


k 


1 

2 


—  *T  — 
+  r,  w 
k 


—  *  T  — 

w  H  w, 
k  k 


(A-4) 


An  optimum  pk  then  can  be  computed  at  each  iteration.  Only  the  choice  of 
w  remains  to  completely  define  the  iterative  process.  It  is  this  choice  of 
Wk  which  distinguishes  the  various  iterative  methods  available.  By  choosing 
Wk  t0  be  the  c°iurtlns  of  identity  matrix  taken  in  cyclic  order,  the  method 
becomes  the  Gauss-Seidel  technique.  This  technique  will  not  be  considered 
in  this  report. 


A.  STEEPEST-DESCENT  METHOD 

If  the  process  is  desired  to  move  in  the  direction  of  ’’steepest 

descent,  "  then  AQ  is  minimized  with  respect  to  the  elements  of  w  .  This 

_  _  _  _  k 

results  in  choosing  w^  =  rk  =  b  -  Hv^.  Substituting  this  choice  into  the 
equation  for  p  yields 


(A-5) 


At  this  point,  p^  is  the  reciprocal  of  the  Rayleigh  quotient  for.  the  residual 
vector.  The  Rayleigh  quotient  has  two  useful  properties: 

•  For  an  arbitrary  vector  r,  the  Rayleigh  quotient 
always  lies  between  the  largest  and  the  smallest 
eigenvalues  of  the  matrix  H 

•  For  a  first-order  approximation  to  an  eigenvector 
of  H,  it  yields  a  second-order  approximation  to  the 
corresponding  eigenvalue 
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These  properties^  imply  that  only  one  calculation  of  pfc  is  necessary  if  the 
initial  choice  of  V()  is  close  to  the  desired  solution  H-l  b.  This  results 


from  the  ability  to  express  the  residual  vector  r 
of  the  normalized  eigenvectors  y  of  H. 


k  as  a  linear  combination 


n 

=  £  atk  vt 

1=  1 


_ ft 

with  y  y  =  1 


(A -6) 


The  matrix  H  can  also  be  expressed  as  the  sum  of  its 
eigenvalues  and  eigenvectors; 


n 


H  =  £  X 
1=  1 


—  — # 
l  Yl  Yl 


for  X.  £  X.  ^  o 
1  J 


(A  -  7) 


This  simple  assumption  is  not  necessary  to  the  conclusion,  and  the  result 
is  still  valid  when  H  is  not  of  this  form. 

Substituting  Equations  A-6  andA-7  into  Equation  A-5  yields 


(A -8) 


Since  orthogonal  eigenvectors  can  be  generated  for  a  Hermitian  matrix, 
Equation  A-8  reduces  to 


Pk  = 


n 


1=1 


atk  a-lk 


n 


^  ^k  alkalk 


(A-9) 
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(A- 16) 


ztk  =  \  rk 


While  the  notation  of  Equations  A-14,  A-15,  and  A-16  is  more 
complicated  than  that  of  Equations  A- 11  and  A- 12,  the  amount  of  computation 
needed  has  been  reduced.  The  reduction  in  computation  is  obtained  by  per¬ 
forming  the  vector-dot-product  operations  x*vk  and  x *  7k  first,  leaving 
only  scalar  operations  to  be  performed.  An  analysis  of  the  amount  of  com¬ 
putation  required  is  presented  in  a  later  section  of  this  report. 

B.  FINITE  PROCEDURE 

By  using  the  method  of  conjugate  gradients,  it  is  possible  to 
choose  pk  and  so  that  the  iteration  process  terminates  in  exactly  n  steps. 
This  is  accomplished  by  allowing  to  be  a  combination  of  the  current  re¬ 
sidual  and  the  previous  vector  w  so  that 

k-  1 


=  ru  +  q 


k- 1  k-1 


(A-  17) 


The  successive  constants  qR_  j  are  chosen  to  make  the  quadratic  foi 


WkHwk-l  =  0 


(A- 18) 


The  constant  pR  is  again  chosen  to  minimize  Q  in  Equation  A-l  so  that 


Pk  = 


—  *  _ 
r  r 
k  rk 

—  *  — 
r,  H  r. 
k  k 


(A-  12) 
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Substituting  Equation  A- 17  into  Equation  A- 18  and  solving  for  the  constant 
q,  l  yields 


q 


k-1 


VkHwk-l 

"k-lH"k-l 


(A- 19) 


The  particular  choices  of  pk  and  qk  ^  just  given  produce 
several  important  effects.  First,  the  vector  product 


rk+i  wk  ■  wk  rk+i 


wk(rkpk 


(A-20) 


vanishes,  which  implies  that  rR+1  and  are  orthogonal.  Second,  the 
vector  product 


— #  — 
rk+l  rk 


=  (rk-  pkwk  H)( 


wk‘  Viwk 


J 


(A-2 1) 


w*  H  w 
k  k-1 


also  vanishes  as  a  result  of  Equation  A- 21  and  the  choice  of  p  .  This  implies 

——  _  k 

that  rk+^  and  rk  are  also  orthogonal.  It  can  be  shown  by  induction  that 

"k+l  *  7k+l  =  7k+l  *  0  (A-22) 

for  l  =  0,  1 ,  .  .  .  ,  k-  1 


That  is,  every  new  rk  computed  is  orthogonal  to  all  the  previously  computed 
residual  vectors.  Similarly,  each  rk  is  orthogonal  to  all  the  previously 
computed  w^  vectors.  It  follows  that  r^  a  0,  since  the  n  orthogonal  vectors 
rQ  through  r^_  ^  span  the  n- space  defined  by  the  matrix  H  of  rank  n.  These 
facts  allow  simplifying  the  computation  of  q  to 

K  “  i. 
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(A-23) 


-*  - 
rk  rk 


k- 1  -*  - 

rk-l  rk-l 


The  computational  cycle  evolving  this  procedure  is 


v,  ,  ,  s  v,  +  p.  w, 
k+1  k  k  k 


(A-  1 0) 


r  =  b  -  H  v, 
k  k 


(A-  1 1) 


Pk  " 


— 

rk  rk 
-*  - 

rkHrk 


(A- 12) 


-*  - 
rk  rk 


k-1  — *  — 

r  r 
k-1  rk-l 


(A-23) 


w,  =  r,  +  q  w 
k  k  4k- 1  k-1 


(A  -  17) 


The  first  cycle  is  computed  using  =  r  . 

Again,  consider  the  special  case  where 


m 


H  =  |I  +  £  a 
1=  1 


-  -* 
l  Xt  Xl 


(A-  13) 
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(A-  10) 

(A- 14) 

(A-24) 

(A-16) 

(A-25) 

(A-26) 


Wk  =  "k  +  Vl  Wk-1  (A- 17) 

For  the  first  iteration,  wQ  =  rQ  and  qk  i  is  not  computed. 

Notation  again  becomes  more  cumbersome,  but  the  computation 
needed  is  much  less  than  that  required  to  form  the  autopower- cross  power 
matrix  and  then  solve  for  the  desired  filter  set.  Computational  savings  are 
obtained  by  carrying  out  the  same  computations  in  a  different  order. 
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APPENDIX  B 

EXACT  INVERSE  EQUATION 
A.  SOLUTION  BY  THE  EXACT  INVERSE  EQUATION 

The  exact  inverse  equation  is  an  analytic  expression  for  the 
inverse  of  the  matrix  H  in  the  equation 

hI  =  b 

The  solution  vector  a  can  be  written  as 

a  =  H'1  b 

when  H  is  nonsingular.  This  report  studies 
Equation  B-l  when  the  H  matrix  is 

m 

H  =  «S+£  a  t  X  *  (B_3) 

l  =  1 


(B-l) 


(B-2) 

numerical  methods  of  solving 


*  Consider  a  matrix  of  the  form  [a  +  x  y^J.  When  both  A  and 

[A  +  x  y  J  are  nonsingular,  the  inverse  of  this  matrix  is 


[a+;7*] 


_1  1  a'1  A-1 

=  a"1  -  A  x  y  A 
— *  - 1  — 
1  +  y  A  x 


(B-4) 


where 


A  =  n  x  n  complex  matrix 


x  =  n  x  1  column  vector 
v  =  1  x  n  row  vector 
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This  equation  is  verified  by  the  following  multiplication: 


A 


-1 


A  x  y  A 

/  —  *  -  1  —  \ 

(i  +  y  a  x ) 


+  *y*] 


(B-5) 


=  A 


-1 


-1 - * 

+  A  x  y 


1 - *  .  i 

A  x  y  A  A 
(l  +7*  A'1  x) 


1 - *  -i - * 

A  x  y  A  x  y 

(>♦7*  a-1  I) 


=  I  +  A 


-1 - * 

x  y 


a-1; 

-*  .-1  - | 
V  -  A  x  ' 

(7 

*  A'1  “'l 
A  x  J 

1  “* 
y 

1 

/  —  *  -  1 
l+y  A 

:) 

1 

=  I  + 


( .  —  *  -  1  —  \  -1 - * 

yi  +  yA  xj  a  xv 


_i  — *  (—*  -i  _i — > 

A  X  y  -  Vy  A  x/  A  xy 


(1,7%-^) 


=  I 


When  the  matrix  H  has  the  form  of  Equation  11-3,  a  recursive 
scheme  can  be  used  to  generate  the  inverse  of  H  directly  from  the  transform 
points.  The  exact  inverse  procedure  finds 

-1 

=  H  (B-6) 


i  S 


m 

1=  1 


—  — * 

\  xl  xi 
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I 


and 


Equation  B-6  is  found  by  first  computing 


=  [*  s + ai  'i  **] 


-i 


(B-7) 


4 


-i 


-l 


a 


6  +ai  xl  S_1  x! 


c-i  -  -*  -1 
S  Xj  Xj  S 


*  [A1  ‘  a2  x2  x2  ] 


-1 


a2‘  ■  a;1 


a 


i+a2x>;lx2 


. - 1  -  -*  -i 
A1  X2X2A1 


a’1  ,-l 

A  *  A 
m  m- 1 


a 


hlJ _ \  A-i  -  -*  -i 

-*  _1  _  I  A  ,  X  x  A 


, 1  +  a  ,  X  A  x 
m- 1  m  m- 1  m 


1  "  "  1 
m-1  m  m  m-1 


-1  .  -1 


H  =  A 


m 


(B-8) 


Each  iteration  is  completed  by  calculating  the  vector 


Cl  =  xi 


Then  the  vector  d  has  the  form 


(B-9) 


(B-  10) 
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The  first  iteration  reviuires  a  knowledge  of  the  matrix  S  and 
the  scalar  £  which  are  stored  in  the  computer  memory  or  calculated  at  the 
beginning  of  each  cycle.  Then,  from  Equation  B-8, 


A 


-1 

1 


a. 


A  + xi  s 


i  - 

X 


-1  -  -*  - 1 
S  x^  x^  S 


(B-ll) 


Each  succeeding  iteration  is 


for  l  =  1,  2,  . .  .  ,  m-  1  (B-  12) 

Since  the  x  's  are  arbitrary  complex  vectors,  this  algorithm  could  be  used 
either  to  stack  a  series  of  autopower-crosspower  matrices  from  successive 
time  gates  or  to  stack  adjacent  frequency  points  (smoothing)  from  the  same 
time  gate. 

B.  HIGH -RESOLUTION  f-k  FILTER  SETS 

The:  high- resolution  technique  uses  an  estimate  of  the  auto¬ 
crosspower  matrix  of  the  form 

H  =  £1  +  z  z  fB-13) 

where 

I  =  n  x  n  identity  matrix 
z  =  n  x  1  column  vector 
z  =  conjugate  transpose  of  z 
£  =  real  constant 


-1 

^  +  1 


=  A, 


-1 


a . 


_ ;V  _ 

1  +  d.  x. 


- * 

c  d 
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The  elements  of  the  vector 
form  vector  x;  i.  e.  , 


z  are  the  normalized  elements  of  the  data  trans- 


!i  =  *i/(*i**) 


f°r  i  =  1,  .  .  .  ,n 


(B-  14) 


Using  the  exact  inverse  equation,  the 


inverse  of  the  H  matrix  is 


H"  1  =  -  ("i  - _ —  2  *  1 


(B- 15) 


Multiplying  bbyH'1  yields  a  filter 


_  ,  - #  — 

a  =  —  _  2  2  b 

*  (*+2*z) 


(B-16) 


Further  simplification  is  obtained  when  I*  z  =  N. 


The  filter  set  becomes 


7  =  r  /  f  \  /  1 


(i)trb) 


- *  _ 

Z  Z  b 


(B- 17) 


Since  the  b  vector  has  one  nonzero  element  in  the  j**1  value, 
Equation  B-17  is  computed  more  simply  as 


*  ■  i(TTT)  ■: 


J  zl 


■  ■  H'  -(tVt) 


^  s  1,2,...,  n 

l  4  j 


(B-  18) 


The  computationai  efficiency  of  this  method  over  all  the  other 
methods  discussed  is  obvious.  The  filter  sets  designed  by  each  technique 

Ufa  WA  mm2  _ A. _ 11  •  1  T 


were  virtually  identical. 
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APPENDIX  C 
NUMERICAL,  RESULTS 

A.  STEEPEST-DESCENT  METHOD 

The  numerical  techniques  discussed  previously  have  been 
applied  to  the  special  problem  of  designing  high- resolution  f-k  spectra 
filter  sets.  In  this  special  situation,  the  techniques  can  be  made  more 
simple  to  achieve  even  greater  numerical  efficiency. 

The  high- resolution  technique  uses  an  estimate  of  the  auto¬ 
crosspower  matrix  of  the  form 

H  =  414  z  z*  (C-l) 


where 


I  =  n  x  n  identity  matrix 

z  =  n  x  1  column  vector 
—  #  _ 
z  =  conjugate  transpose  of  z 

4  =  real  constant 

The  elements  of  the  vector  z  are  the  normalised  elements  of  the  data  trans¬ 
form  vector  x;  i.  e.  , 


z. 

x 


In  this  situation,  the  residual  vector  r  is 

k 


for  i  =  1,  .  .  .  ,  n  (C-2) 


=  b 


•  [$I  +  Z  z*~\  \ 


(C-3) 
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o 


The  iteration  factor  pk  is  not  affected  by  a  first-order  change 
in  the  residual  vector  r^.  Further  observations  about  the  iteration  factor 
can  be  made  as  a  result  of  the  special  form  of  the  H  matrix  being  considered. 

It  is  desired  to  compete  the  iteration  factor 


_*  _ 
r  H  r 
k  k 


(C-4) 


The  related  eigenvalue  problem  is 


Hyi  *  \y. t 


for  y*  y  =  1 


(C-5) 


where 


l  =  1. 2,  .  .  .  ,  n 


=  eigenvalue  of  H 
=  the  related  eigenvector 


- * 


When  H  is  (4  I  +  x  x  ),  it  has  only  one  eigenvalue  X  , 

n 


—  *  — 

X  =  £  +  x  x 
n 


(C-6) 


and  all  the  n- 1  other  eigenvalues  are  £. 


Suppose  an  initial  vector  is  chosen  so  that 


r.  =  a  y  +  s 
1  n 


(C-7) 
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where  ^  is  the  eigenvector  corresponding  to  the  nonzero  eigenvalue  X  and 

e  is  a  vector  representing  the  difference  between  7  and  a  y  .  Then 

1  n  ’ 


p  .  [ayn%e1KH 

[ay*  +  r*]H[a7n +  F] 


(C-8) 


l  —  *  —  r-*  -  -*  _T  _*  _ 

-  a  ynyn+a[y  Y  l+  £  £ 

2  —  '>c  —  r  —  #  —  *  -I  vt, 

a  yn  Hyn+aLyn  He  +  c  H  YnJ  +  6  H 


By  replacing!^  with  Xr  yn  in  the  denominator  of  Equation  C-8  and  recog¬ 
nizing  that  yn  y  =  1  by  definition,  it  is  possible  to  write 


2  /— *  —  — 
a  +  ai  y  e  +  €  y  )  +  e  e 

Pl  =  — - _ jhi _ 

\i(.a  +  aXyn  ®  +  ®  yn)J  +  e  H 


(C-9) 


If  the  difference  vector  e  is  small  compared  to  the  eigenvector 
yn’  then  Equation  C-9  can  be  approximated  as 


pi  =T 


(C-10) 


=  4  +  x' 


then. 


=  (-^) 

'  i  +  x  x  / 


(C-  1 1) 
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When  the  H  matrix  has  the  form  in  Equation  C-l, 


P 


1 


1 

£  +  n 


(C-  12) 


where  n  is  the  number  of  elements  in  the  vector  z.  The  value  of  p^  is  observed 
in  Figure  C-  1 . 

The  objective  then  is  to  find  an  initial  vector  v^  which  will  re¬ 
sult  in  a  small  e  component  in  the  first  residual  vector.  If  this  is  found,  the 
iteration  factor  need  not  be  computed  at  all  but  can  be  taken  to  be  the  value 
of  Equation  C-12. 

High- resolution  f-k  spectra  are  generated  from  a  filter  set  a 
designed  using  an  H  matrix  in  the  form  of  Equation  C-l  and 

o  =  [l,  0,  0,  .  .  .  ,  oj  (C- 13) 


The  unit  entry  in  b  is  placed  in  the  position  corresponding  to  the  desired 
reference  sensor.  All  other  elements  of  the  b  vector  are  0.  The  equation 
set  being  solved  is  then 

r  ^ 


(C  -  14) 


0 


v. 
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Figure  C  -  1 .  Variation  of  from  Iteration  to  Iteration  when  the 
_ Residual  Is  Small  (Real  Data;  £  =  0.001) _ 
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Experimental  results  show  that  =  b  is  a  good  choice  for 

the  initial  ’’guess"  at  a.  Using  this  value  for  v^,  Equation  C- 12  to  solve 

Equation  C- 14,  and  the  steepest-descent  method,  an  accurate  answer  often 

was  reached  in  one  iteration;  seldom  were  more  than  three  iterations  re¬ 
ft 

quir^d.  Bodewig  has  mentioned  that  0.9  p^  sometimes  produces  better 
convergence  than  p^.  On  real  data,  0.  9  p^  is  not  so  good  a  choice  as  p  . 
The  convergence  rates  using  different  iteration  factors  can  be  observed  in 
the  Figure  C-2  plots  for  each  iteration  of  error  vs  iteration  constants. 

When  synthetic  data  were  used,  the  0.9  p^  iteration  factor  converged  faster 
than  Pj  (Figure  C-3). 

If  a  good  approximation  to  the  filter  weights  is  not  available, 
very  slow  convergence  can  result.  This  is  demonstrated  in  Figure  C-4  for 
synthetic  data  with  an  initial  vector  v^  value  of  (1  +  jl)  for  every  entry 
where  j  =  /T. 

There  is  an  indication  that  solving  the  degenerate  problem 
caused  the  rapid  convergence  experienced  using  a  steepest-descent  design 
procedure.  An  initially  good  "guess"  at  the  filter  weights  is  not  avail¬ 
able  when  the  H  matrix  contains  the  sum  of  many  vector  products,  implying 
that,  in  general,  relatively  slov  convergence  can  be  expected  in  the  more 
complicated  problems. 


B.  FINITE  PROCEDURE 

Results  also  have  been  good  using  the  finite  procedure.  The 
chief  advantage  of  this  method  —  that  it  is  not  necessary  to  start  from  a 
first  approximation  to  the  filter  weights  as  was  necessary  when  using  the 
steepest-descent  method  —  is  shown  in  Figure  C-5. 


Bodewig,  E.  ,  1956:  Matrix  Calculus,  Interscience  Publishers ,  Inc., 
New  York. 
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ing  a  p1  Which  Is  Slightly  J 
>est  descent;  real  data;  £  = 


CHANNELS  2  THROUGH  4 


CHANNELS  5  THROUGH  8 


Pi  =°-9popt 
Pi  =  Popt 


CHANNELS  9  THROUGH  12 
CHANNELS  13  THROUGH  16 


-0.  335 


-0.460 


-0-0.  503 


-0.552 
-0.  570 


NUMBER  OF  ITERATIONS 


Figure  C -3.  Approximate  Convergence  of  Filter  Weights  When  Syntheti 
Data  Are  Used.  (Plane-wave  incident  on  16-channel  square 
array;  channel  1  is  reference;  £  =  0.  1) 


TO  0.  943 


Figure  C-4.  Slow  Convergence  of  Filter  Weights  When  All  Initial  Filter 
Weights  Are  1  (Synthetic  data;  £  =  0.  1) 
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CHANNELS  2  THROUGH  4 
CHANNELS  5  THROUGH  8 
CHANNELS  9  THROUGH  12 
CHANNELS  13  THROUGH  16 


CHANNEL  1  IS  REFERENCE  (NOT  SHOWN) 

- L - 1 _ I _ I _ I 

3  4  5  6  7 

NUMBER  OF  ITERATIONS 


Figure  C-5.  Approximate  Convergence  of  Filter  Weights  When  All  Weights 
Are  Initially  Set  Equal  to  1.  [Finite  procedure;  synthetic  data 
£  =  0.  1;  (£  I  +  x  x*)  a  =  ‘Bl 


Using  real  data,  convergence  to  essentially  the  exact  solution 
was  obtained  in  two  iterations;  every  value  of  the  initial  filters  was  1  (Fig- 
ure  C-6)-  Of  course,  when  the  real  part  of  the  reference  channel  was  set 
to  1  and  all  other  weights  set  equal  to  0,  much  better  convergence  occurred 

(Figure  C-7).  In  fact,  the  filters  were  correct  to  four  significant  figures 
after  the  first  iteration. 

Since  the  imaginary  part  of  p^  was  very  small  during  the 
first  few  iterations  (Table  C-l),  the  program  was  rewritten  to  set  the 
imaginary  part  to  0,  tnereby  reducing  the  number  of  computations.  The 
convergence  was  completely  unaffected. 

Finally,  p^  =  0.  9  p^  was  used  to  measure  the  effect  of  varying 

Pk'  FUter  wei&hts  diverged  rapidly  (Figure  C-8);  the  strong  susceptibility 

to  changes  in  p  indicated  that  p  necessarily  must  be  recalculated  at  each 

iteration,  especially  in  view  of  the  large  and  apparently  unpredictable 

(without  recalculating)  variations  in  p  (Table  C-l). 

k 

Although  the  number  of  computations  involved  is  almost 
twice  that  required  by  the  steepest-des  cent  method,  the  finite  procedure  is 
better  than  any  direct  method  other  than  the  exact  inverse  equation.  The 
rapid  convergence  of  this  method  is  certainly  the  result  of  the  degenerate 
problem  being  considered. 

C.  SPECTRA 

Wavenumber  spectra  obtained  from  the  iteration  schemes 
presented  compare  quite  favorably  with  those  obtained  from  current  high- 
resolution  techniques.  Spectra  were  obtained  from  real  data  (C3  subarray, 
minimum  velocity  10.  0  km/sec,  f  =  0.  £5  Hz)  using  filter  weights  designed 
by  the  exact  inverse  equation  and  by  the  finite  procedure.  In  both  cases, 
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Table  C-l 


VARIATION  OF  p  IN  FINITE  PROCEDURE  (REAL  DATA) 


Iteration 

p  (real  part) 

p  (imaginary  part) 

1 

-0. 0400 

0.  ; 

2 

-0. 0400 

2.50  x  10"10 

3 

-0. 0400 

7.54  x  10‘6 

4 

-0. 0278 

0. 0346 

5 

0. 0104 

Note:  Initially,  the  real  part  of  the  reference  channel  is  1  and  all  othe: 
channels  are  0. 


Iteration 

p  (imaginary  part) 

1 

-0. 0400 

0. 

2 

-999. 

-2.76  x  10  8 

3 

-0. 0399 

-9.63  x  10"10 

4 

-0. 0417 

-4.31 x  10-8 

5 

-0. 0244 

-4.  78  x  10-4 

6 

-0. 0109 

-1.83  x  10'6 

■he  filter  weigh,  of  channel  8  has  a  large  because  ^  was 

punched  incorrectly  on  the  data  card.  However,  the  spectra  obtained  are 
virtually  identical  with  those  obtained  from  current  techniques.  Since  the 
error  in  channel  8  did  no,  significantly  affect  the  spectra,  wavenumber 
spectra  are  seen  to  be  relatively  insensitive  to  errors  in  the  individual 
i  ter  weights  (unless  the  error  happens  to  be  in  the  reference  channel). 

Filter  weights  designed  from  the  steepest-descen,  method 
were  virtually  identical  to  those  obtained  from  the  exact  inverse  equation; 

ere  ore,  the  steepest-descent  method  also  will  yield  good  spectra. 

D.  SMOOTHING 

Fairly  good  results  were  obtained  when  data  were  smoothed 
over  a  range  of  frequencies.  The  data  used  were  derived  from  the  data 
for  the  firs,  eight  channels  of  the  real  data  from  the  C3  LASA  subarray 
The  transform  for  a  given  channel  was  altered  to  give  four  new  transforms 
with  an  average  equal  to  that  of  the  original  transform.  The  deviation  of 
any  given  component  from  that  of  the  transform  from  which  i,  was  derived 
was  1  to  2  percent.  If  these  deviations  were  truly  random,  the  smoothed 
i  ter  weights  would  be  expected  to  be  equal  to  those  for  the  original  data. 

Since  the  deviations  were  no,  entirely  random,  however,  they  were  actually 
somewhat  different. 

With  all  initial  filter  weights  set  to  0  except  for  the  real  part 
o  the  reference  channel  which  was  set  to  1,  the  steepest-descent  method 
showed  a  large  initial  correction  followed  by  very  slow  convergence.  This 
method,  therefore,  seems  to  be  somewhat  ineffective. 

The  finite  procedure,  on  the  other  hand,  produced  convergence 
to  a  nonvarying  set  of  filter  weights  after  three  iterations  and  was  reasonably 
correct  (two  significant  figures)  after  two.  These  results  occurred  only  when 
the  initial  filter  se,  was  a  firs,  approximation  to  the  correct  se,  of  weights. 
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Whe  a  bad  first  approximation  was  made  (e.g.  ,  letting  the  imaginary  par, 

the  reference  channel  be  se,  to  .),  very  slow  convergence  followed  a 
large  initial  correction. 

ntatrice.  ^  ““  ^  Pr°“d"”  Pr°du“s  «°°<1  -suits  on  full  rank 

3  reaS°nable  approximation  to  the  filter  weights  is  available  The 

rrscenr ,hod  does  not  — 

1  the  filter  weights  is  available. 

.  ^  Thfc  necessity  of  a  good  first  approximation  is  further  illus¬ 

trated  when  data  from  the  10-channel  CPO 
,  ,  .  1161  CPO  arraY  were  smoothed  over  a  range 

requencies;  (only  nine  channels  were  used)  The  h  t 
-  .  used).  The  data  sample  consisted 

Of  nets  1.00  Hr  for  ,6  October  1,64,  with  41  Pourier  transform  taken 
•  95  to  1. 05  Hz  at  frequency  increments  of  0.  0025  Hz,  All  trans 
orms  were  equally  weighted.  Neither  method  could  produce  convergence 
when  the  real  part  of  the  reference  channel  was  1, 0  or  0.  1.  Table  C-2 

show,  the  filter  weights  actually  obtained  using  other  methods.  I,  is  readily 

apparent  that  the  initial  choices  HiH 

dices  did  not  approximate  these  filter  weights 

an  also  that  no  general  way  to  make  a  good  approximation  seems  apparent. 

should  be  noted  that  very  poor  resolution  was  obtained  with  the  filter 

weights  of  Table  C-2,  indicating  that  the  data  were  very  poor. 

Effective  iterative  methods  for  smoothing  are  thus  limited 

o  case,  in  which  the  transforms  at  the  various  frequencies  under  consider- 

ion  are  approximately  equal,  allowing  all  initial  filter  weights  to  be 
approximated  as  0  (excpnt-  th* 

that  th.  r  -1  "nel  Which  U  lt  appears 

that  the  finite  procedure  is  preferable  to  rh  4 

e  steepest -descent  method 

when  computing  spectra  for  smoothed  crosspower  matrices. 
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Table  C-2 

FILTER  WEIGHTS  FOR  DATA  SMOOTHED  OVER  41  FREQUENCIES 

FROM  0.95  TO  1.05  Hz 


Channel 

Filter 

Weights 

1 

0.  22481  x  10‘10 

-0. 71687  x  10'7 

2 

-0.41205  x  10"2 

0. 16657  x  10'1 

3 

0.  29080  x  10"2 

0.  5448  x  10"2 

4 

-0. 12049  x  10'2 

0.  16102  x  10"2 

5 

-0. 10966  x  10 

-0. 13375  x  10"3 

6 

-0.78141  x  10‘3 

0.84145  x  10"4 

7 

0. 14035  x  10_1 

0.40036  x  10"2 

8 

-0. 96124  x  10"2 

-0. 75366  x  10"2 

9 

0.  22133  x  10“2 

-0.96892  x  10"2 
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